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Abstract 

We study the polydisperse Baxter model of sticky hard spheres (SHS) in the modified 
Mean Spherical Approximation (mMSA). This closure is known to be the zero-order ap- 
proximation (CO) of the Percus-Yevick (PY) closure in a density expansion. The simplicity 
of the closure allows a full analytical study of the model. In particular we study stabil- 
ity boundaries, the percolation threshold, and the gas-liquid coexistence curves. Various 
possible sub-cases of the model are treated in details. Although the detailed behavior 
depends upon the particularly chosen case, we find that, in general, polydispersity in- 
hibits instabilities, increases the extent of the non percolating phase, and diminishes the 
size of the gas-liquid coexistence region. We also consider the first-order improvement of 
the mMSA (CO) closure (CI) and compare the percolation and gas-liquid boundaries for 
the one-component system with recent Monte Carlo simulations. Our results provide a 
qualitative understanding of the effect of polydispersity on SHS models and are expected 
to shed new light on the applicability of SHS models for colloidal mixtures. 
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I Introduction 



In sterically stabilized colloidal mixtures, particles are coated with polymer brushes to prevent 
irreversible flocculation due to van der Waals attraction [1,. If the solvent is a moderate one, 
a lowering of the temperature yields very strong attraction with a range much less than the 
typical colloidal size. In microemulsions of polydispersed spherical water droplets each coated 
by a monolayer of sodium di-2-ethylhexylsulfosuccinate dispersed in a continuum of oil, the 
droplets interact with each other via a hard core plus a short range attractive potential, the 
strength of which increases with temperature |2j- For these systems, a very useful theoretical 
model is the sticky hard sphere (SHS) model proposed by Baxter [3] long time ago for atomic 
liquids. In the original Baxter solution [2 E] the one-component Or stein- Zernike (OZ) integral 
equation was analytically solved within the Percus-Yevick (PY) approximation. Successive 
extension to mixtures pij, however, proved to be a formidable task in view of the fact that a 
large (infinite 1 ) number of coupled quadratic equations ought to be solved numerically in order 
to have a complete understanding of both thermodynamics and structure of the model. This is 
the reason why, to the best of our knowledge, only binary mixtures have been explicitly discussed 
so far in this framework Moreover it has been proven by Stell jHJ that sticky spheres of 
equal diameter in the Baxter limit are not thermodynamically stable and size polydispersity 
can be expected to restore thermodynamic stability. 

Motivated by this scenario, it was recently proposed 7\ a simpler approximation (mMSA 
closure) having the advantage that also the multicomponent case could be worked out analyt- 
ically [HI Ej. Further analysis and comparison with both Monte Carlo (MC) and PY results 
[3 EE E] i n the one-component case, have shown that the mMSA closure for Baxter model is 
a reliable one up to experimentally significant densities. The price to pay for this simplification 
is that only the energy equation of state gives rise to a critical behavior, the other two routes 
yielding either a non-critical behavior (compressibility), or a diverging equation of state (virial). 

In this work we pursue this investigation by studying the multicomponent version of the 
model proposed in Ref. [Zj, and analyzing various consequences. We first solve the multi- 
component version of Baxter model within the mMSA closure, and show that the solution is 
equivalent to the one derived in Ref. [8 for a companion SHS model. The solution, derived in 
terms of an auxiliary function called Baxter factor correlation, turns out to be formally similar 
to that derived with the PY closure. However, and this is the crux of the matter, the matrix 
function representing the stickiness parameters is unconstrained, unlike the PY counterpart. 
In order to make further progress and derive the multicomponent energy equation of state, a 
further assumption is necessary on the matrix representing the stickiness parameters. As dis- 
cussed previously (see Ref. jH] for details) a remarkable simplification occurs when the general 
element of this matrix has the form of a sum of dyads (i.e. it is dyadic). In these cases the 
necessary matrix inversion can be carried out analytically and all measurable quantities can 
then be computed. Physically, this reduction to a dyadic form amounts to assume a relation 
among polydispersity in size and polydispersity in stickiness, that is on the adhesion forces. In 
addition to the two cases proposed in Ref. [H] (denoted as Case I and II in the following) and 

1 Strictly speaking we should distinguish between discrete polydispersity (multicomponent mixture with a 
large number of components p f=a 10 2 -j- 10 3 ) and a continuous polydispersity corresponding to p — > oo with a 
continuous distribution of sizes or other properties. This distinction will be specified in more details in Section 
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that proposed in Ref. [JJ] (Case IV), we shall consider two further cases. The first one (Case 
III) is a physically motivated variant of Case I, whereas the second one (Case V) has its main 
justification in the simplifying features occurring when one attempts to go beyond the mMSA 
closure with a density perturbative approach (to first order this will be called CI, as in Ref. 
[7], for reasons which will become apparent in the rest of the paper). 

The main results of our analysis are the following. We derive the instability curves in three 
of the considered cases (Case I-III) within the mMSA approximation and analyze the effect of 
polydispersity in some detail. In order to test the reliability of the mMSA approximation, we 
also consider the first-order correction (CI) in the one-component case and compare with the 
PY result. 

Next we consider the effect of polydispersity on the percolation threshold. This is an 
interesting phenomenon on its own right and has attracted considerable attention recently 
[T ^ lT H lT5 | IT! H lllj. being a paradigmatic example of flocculation instability. In particular, recent 
Monte Carlo simulations pTHITT] on monodisperse (one-component) spheres with sticky adhesion 
have clearly tested the performance of analytical calculations based on the PY approximation 
[Til I15j . We then study the percolation transition as a function of polydispersity in all above 
mentioned cases within mMSA. Again we can discriminate the effect of polydispersity on the 
percolation line, and also compare it with the first-order correction CI, the PY approximation 
and MC simulations in the one-component case. 

Next we consider phase equilibrium. A major obstacle to the analysis of phase transition 
in polydisperse systems is posed by the fact that, in principle, one has to deal with a large 
(infinite) number of integral non-linear equations corresponding to the coexistence conditions 
among various phases. In this model, however, as it also occurs in other simpler models such as 
hard spheres (HS) ^B] , van der Waals fluids ^JJ and in more complex cases such as factorizable 
hard-sphere Yukawa potentials [1111111, the task can be carried out in full detail in view of the 
fact that the (excess) free energy depends upon only a finite number of moments of the size 
distribution function. In the particular case of two-phase coexistence, we derive the cloud and 
shadow curves of all Cases in the mMSA approximation. We compare the results with those 
derived earlier for a polydisperse van der Waals fluid 17, and discuss analogies and differences 
in this respect. Finally we compare the results of the mMSA one-component case with the 
first-order correction, the PY approximation, and the results of MC simulations. 

The plan of the paper is as follows. In Section [hi] we define the multicomponent SHS model, 
give the solution for Baxter factor correlation function in the mMSA (CO) approximation, and 
define the various Cases of polydispersion models taken under exam; In Section Um we give the 
solution for Baxter factor correlation function in the CI approximation and show how Case V is 
particularly suitable to study the polydisperse system analytically; in Section [TV] we analytically 
derive the instability boundaries; in Section [3 we find analytically the percolation thresholds; 
In Section IV1I we derive numerically the two phase coexistence curves; In Section IV11I we lay 
down our conclusions and further developments. 
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II Baxter model and modified MSA solution 



In Baxter model of sticky hard spheres (SHS1), one starts adding to the hard sphere (HS) 
potential a square-well tail with 



Mr) = -k B T\n ( 1 ) , a i3 < r < R l3 , (2.1) 

V U ^3 ^3 / 

where a i3 = (<jj + &j)/2 (crj being the HS diameter of species i), R i3 — a i3 denotes the well 
width, kg is Boltzmann constant, T the temperature, and the dimensionless parameter r^ 1 > 
measures the strength of surface adhesiveness or 'stickiness' between particles of species i and 
j (jij is also an unspecified increasing function of T). The sticky limit corresponds to taking 
{Rij} -> {rr, ; j. 

The Baxter form of the Ornstein-Zernike (OZ) integral equations for this model admits a 
very simple analytic solution if one uses the following modified Mean Spherical Approximation 
(mMSA) 

c ij ( r ) = fij( r ) for r > Vij , (2-2) 

where (r) and /^(r) = exp [—(3<f>ij(r)] — 1 are the direct correlation function and the Mayer 
function, respectively [f3 = (/c^T) -1 ]. This can be easily inferred by using the formalism 
introduced in Ref. [7j. As pointed out in that reference, the mMSA closure can be reckoned as 
a zero-order approximation in a perturbative expansion, and hence it will also be denoted as 
CO henceforth. In terms of Baxter factor correlation functions % (r), its extension to mixtures 
reads 

( \ _ / \ a i( r ~ a ij? + ( b i + a i a ij)i r ~ + K ij , L ij = iPi ~ °j)/2 < r < (Tij , , , 

qi3[T> ~\0, elsewhere , 1 6} 

1 36^ 120 h [I \<Ti (0A s 

ai = A + ~& A~ ' bi= \-- ai )-> lL[} 

p p 



A A 2 A ' \A 72 

7T P 

in = Pt^i , Ci=-J] p m a m K im , A = 1 - £ 3 , (2.5) 



^-f 6 

«=1 m=l 

with p being the number of components, pi the number density of species i, and 

j^faMSA) = _^_ a% s K% p 6) 

We remark that although Eqs. ()2.3)) - ([2.5|) are formally identical to their PY counterpart, this 
result is in fact simpler in such they differ in the quantity K i3 which in the PY approximation 
reads [H] 

where yij Y \<Jij) is the contact value of the PY cavity function. In general, the parameters \j 
can be determined only numerically by solving a set of p(p+ l)/2 coupled quadratic equations 
|2*m |3] , and this makes the multicomponent PY solution of limited interest from the practical 
viewpoint. In particular a global analysis of the phase diagram proves to be a formidable task 
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within the PY approximation [3]. On the other hand, in view of the simplicity of Eq. ()2.6|) 
with respect to its PY counterpart Eq. (|2.7|) . this is indeed possible within the mMSA (CO) 
approximation. The above results is, moreover, fully equivalent to a parallel but different sticky 
HS model (SHS3) studied by us in previous work jHUH]. Hence, as discussed in those references, 
this analysis can be pursued analytically provided that Kij has a dyadic form. To this aim, we 
consider polydisperse fluids with HS diameters distributed according to a Schulz distribution 2 . 

As regards stickiness, we choose to keep it either constant or related to the particle size. 
There are two main reasons for this. First, one expects the adhesion forces to depend upon the 
area of the contact surface between two particles (see Fig. |TJ), and hence on their sizes. Second 
and more practical reason, is that this is a simple way of obtaining the required factorization. 
As the stickiness-size relation is not clearly understood, we consider five different possibilities, 
denoted as Case I-V henceforth. The three simplest choices are 



1 (a) 



2 



Tij t a% ' 
1 _ 1 a.aj 



j^(mMSA) 



J_ = l <^> 
Tij T (7?- ' 



(mMSA) 



^■{mMSA) 



Case I 12r 
1 



' <^> 2 , (2.f 



Case II 12r 
1 



Case III 12r 



> (2-9) 
<a 2 ) . (2.10) 



where (a) is the average HS diameter ((F) = ^XjFj, here Xi = pi/p is the molar fraction of 
species i with p = £\ pi the total number density) and r is assumed to depend only on the 
temperature, while the remaining factor in t^ 1 is a measure of stickiness strength and is related 
to the particle sizes. The physical interpretation of these choices is the following. In Case I the 
stickiness is assumed to be proportional to the surface contact area of two colloidal particles 
having average size (a), whereas in Case II the adhesion of each particle is linearly related to 
its size. Case III, finally, is a variant of Case I where one considers an average stickiness rather 
than the stickiness of an average particle. 

In all these cases the j{^ MSA > matrix can be factorized as 

K ^ MSA) = Y l Y j , (2.11) 

with Yi having dimensions of length (Yi = (Vl2r) 1 (a), Yi = (Vl2r) 1 <7j, and Yi = 

(\/12t) (a 2 ) 1 ^ 2 in Case I, II, and III, respectively). Note that Case I and II have already 
been exploited by us in previous work jS]. 

We also consider a case similar to that proposed by Tutschka and Kahl ^2] (henceforth 
denoted as Case IV) 

11 , 
— = - , (2.12) 

Tij r 

In this case the j^( mM5A ) matrix can be written as a sum of three factorized terms (as it can 
be immediately inferred by expanding the square a 2 - = (o"j + (Tj) 2 /4) and has the interesting 

2 Here, for simplicity, we disregard possible complicancies arising from the fact that unphysically large par- 
ticles are included in this analysis. These were discussed in Ref. |18| . 
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physical interpretation of being proportional to the area of the actual contact surface 47rcr?- 
between particles of species % and j. Finally, and for reasons related to the CI approximation 
that will be further elaborated below, we consider Case V defined by the linear (rather than 
quadratic) dependence 

1 1 (a) 



;2.i3) 



in this case the j^^- MSA ) parameters can be written as a sum of two factorized terms. 



Ill The CI approximation 

It was recently argued [7j in the one-component case, that the mMSA (CO) approximation can 
be improved by including the next order term in the density expansion of the direct correlation 
function. Its extention to multicomponent mixtures reads 

Cij{r) = fq{r)[l + $> m 7^(r)] r > ay , (3.1) 

where 

7$,(r) = J f im (\r-r'\)f mj (r')dr' 

= — / dssf im (s) / dttf mj (t) . (3.2) 

r J J\r-s\ 

is the first-order coefficient in the density expansion of the partial indirect correlation functions 
~fij(r). As discussed in Ref. [7], if we retain in the PY closure only the terms corresponding to 
the zero and first-order expansion in density we recover the CI approximation (j3.1|) . It turns 
out that Baxter factor correlation function can still be cast in the form, Eqs. ()2.3|) - (j2.5|) but 
the K{j parameters have the form 

= K^y^iati) , (3.3) 
where the partial cavity functions at contact for this closure are 

i/r^^+EwSK-) , (3.4) 
in 

Using in Eq. ()3.2|) fij(r) = —Q{oij — r) + 8(r — ay)crij/(12r^), we find after some algebra 
the following result 

r 2 



2tt I o im 



2 

\(a 2 - J 2 \ 4- °" mj 



mj J 



+ 



9 a 2 • 1 

n 0i 3 mi i q V^mi mi) ' 

O JL — / j _ 

~^y (J mj ~ °'y)( £7 'mi — ^mi) T^^iji^mi ~ ^mi) ~ 



Because of the presence of the factor l/ay in Eq. ([3.5)1 . K-?^ cannot be expressed as a sum 
of factorized terms if we use any of the Cases I, II, or III. Case IV, on the other hand, would 
be tractable, but it would yield K^ 1 ' as a sum of 14 factorized terms (proportional to a^aj 1 
with n,m = 0, 1, 2, 3 except n = m = 0, 3) which is unmanageable in practice. In Case V, on 
the other hand, a great simplification occurs and we find 



where 



K ij 1] = k o + (o-j + 0j)ki + a i a j k 2 , (3.6) 



1 (a) 3 (a 2 ) 1 

k ° = ^sre^T^' (3 - 7) 

11 / 1 (a) 4 1 1 (a) 2 (a 2 ) 1 1 . x 1\ 

iWi (39) 



576 ((T 3 )r 3 24 (a 3 ) r 2 8 (a 3 ) r, 

where rj = £3 is the packing fraction. The expression (J3.6)) is slightly more complicated than 
the x\^ MSA ' treated with Case IV, because of the k term. This noteworthy feature is the 
main justification for the particular form of Case V. 



IV Phase instabilities 

Our first task is the analysis of the phase instabilities for the polydisperse system only in the 
mMSA using Cases I, II, and III. 

The next level of approximation (CI) is considerably more laborious (since the calculations 
for the CI approximation even in the simple case of Case V requires determinants of n-dyadic 
matrices with n > 4) and we shall limit ourselves to the one-component case for simplicity. 



IV. 1 mMSA approximation for the discrete polydisperse system 

For p-component mixtures, one can define the following generalization of the Bhatia-Thornton 
concentration-concentration structure factor 

Scc(k) / (]>-) = l S ^)l E X'^ > 

where |S(/c)| denotes the determinant of the matrix S(k) whose elements are the Ashcroft- 
Langreth partial structure factors |23j. Furthermore, the S^ 1 (k) functions are the elements of 
the inverse of S(k), which can be expressed as 

(*) = - M 1/2 c l3 (k) = J2 Qmi (~k) Q mj (k) , (4.2) 

m 

with Cij (k) three-dimensional Fourier transform of (r), Qij (k) = 5^ — 2-K(p i pj) 1 ^ 2 q i j (k) , 
and %j (k) being the uni-dimensional Fourier transform of (k is the magnitude of the 

exchanged wave vector, 5y the Kronecker delta). 
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Phase instability corresponds to the divergence of the long wavelength limit Sqc (k = 0), 
which is related to the concentration fluctuations. Taking into account the relations 



j2(^) 1/2 sr J l (o) = j2^ = (pkBTK T r i = (?£p) , 

in A \ I / T 



(4.3) 



|S(0)| = |I — c 



Q(o) 



(4.4) 



[where Kt is the isothermal compressibility, I the unit matrix of order p, and C has elements 
(piPjY^Cij (k)], Sqc (k = 0) can be re-expressed as 



Sec (0) 

Elm X m 



1 



Q(0) ( P k B TK T ) 



(4.5) 



For a one-component system the divergence of Kt signals mechanical instability, associated 
with a gas-liquid phase transition or condensation. However, a multi-component fluid usually 
becomes unstable while Kt remains finite and different from zero. In this case, it is the 



vanishing of 



Q(o) 



which causes the divergence of 5cc (0) and produces a phase instability 

(221 122] ■ Indeed if one tries to calculate the locus of points in the phase diagram (r, rj) where 
Yli x i a 'i = 0; using Cases I, II, or III, discovers that such curves disappear (the quadratic 
equations in r have a negative discriminant) as soon as we switch on the size polydispersity 
letting {a 2 ) ^ (o - ) 2 - We remark that the exact nature of this instability requires a more involved 
analysis and it will be deferred to a future work. 



The computation of 



Q(o) 



which usually becomes a formidable task with increasing p, is 



rather simple for the mMSA solution of Baxter model when K^ is factorized as in Eq. (j2.11|) . 
In fact, Q(k) becomes a n-dyadic (or Jacobi) matrix 



[I, J 



(4.6) 



with the remarkable property that its determinant, which is of order p, turns out to be equal 
to a determinant of order n (-C p for polydisperse fluids) [Hj. The necessary expressions are 
reported in Appendix IVHI 

For factorized K^s, one finds 



Qij (o) 



r 7T 

°ij + g [PiPj 



,1/2 



1 3 

A J 



with 



A 



1 



p (a m Y n ) 



, (4-7) 



(4.8) 



((•••) denotes a compositional average, i.e. (FG) = ^2 i x i FiG i ). Note that £ mi0 = ^ 



We emphasize that the decomposition of Eq. ()4.7j) into Af and Bp is not unique. However, 
Qij (0) of Case I and III is 3-dyadic (i.e. it contains n = 3 dyadic terms), while Q^ (0) of Case 
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II is simply 2-dyadic. As a consequence, one has to calculate at most a determinant of order 3. 
The general result for all three cases is 



Q(°) = ^2 K 1 + 2 ^ (1 - 12 6,2) + 36£ 2 2 J 



Physically admissible states must satisfy the inequality 



Q(o) 



> o 



(4.9) 



and the stability 



boundary is reached when 



Q(o) 



0, which yields 



(*)(o*)\ 2 3 V 2 



v(i-v) 

1 + 2/7 

/ ax V 

(a 3 ) 



(a 3 ) 2 1 + 2 V 



Case I , 
Case II , 
Case III . 



(4.10) 



If the HS diameters follow a Schulz distribution, then the stability boundary of Cases I and III 
can be expressed as 



1 



n 



3r) 



M X M 2 
1 Mi 

Mo' 



Mf 



1 + 277 
3i] 



M| 1 + 2?7 



Case I , 
Case III 



(4.11) 



where Mj = 1 +js 2 with s = [(a 2 ) - (a) 2 ] / (a) measuring the degree of size polydispersity. 

The fluid is stable at 'temperatures' r higher than those given by the previous equations 
(since |Q(0)| > 0). Let us now compare two mixtures with the same packing fraction rj but 
different polydispersity degree s. As depicted in Fig. [2] at small r\ values, increasing s at fixed 
7] lowers the stability curve of Case I and III. As shown by the left branch of the curve (the 
opposite trend on the right hand side of the figure is not acceptable, since the mMSA closure 
can be a reasonable approximation only in the low density regime) the onset of instability 
occurs at lower r. As expected, polydispersity renders the mixture more stable with respect to 
concentration fluctuations. Quite surprisingly, on the other hand, the stability boundary does 
not depend on s at fixed r] in Case II, and all mixtures with different polydispersity have the 
same stability boundary as the one-component case (s = 0). 



IV. 2 CI approximation for the one component system 



As remarked, the CI approximation yields rather more complex expressions and here we restrict 
to the one-component case. Yet, this example provides a flavor of how this approximation would 



Q(o) 



0. 



work in the multicomponent case and could be compared with the result given by 

For the one-component system phase instability coincides with the divergence of Kt- As from 
Eg. (|OD 



(pk B TK T y 



V 



T ' 







(4.12) 
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where [see Eqs. ()3.4|) and ([3.5p ] 

y {C1 \a) = l + y 1 {r)r ] , (4.13) 

with 

The curve for the onset of mechanical instability is shown in Fig. El and compared with the 
PY one 

10 -9/(1-,) + 14, 

12(1 + 2ri) K J 

One clearly sees that the CI stability boundary lowers and shifts to the left in agreement with 
the PY result. 



V Percolation threshold 

In view of the simplicity of the mMSA (CO) solution, one might expect that other quantities, 
besides those discussed so far, can be computed analytically. We now show that this is indeed 
the case. The problem we address in this section is continuum percolation. This problem is far 
from being new j2Hj- However new activity along this line has been stirred by recent and precise 
Monte Carlo results for the one-component case [TUt ITT] , and it is then rather interesting to 
consider its multicomponent extension. For the sake of completeness we now recall the basic 
necessary formalism ^Hl El 02] ■ 

In the sticky limit the partial Boltzmann factors read 

K°- 

13 



(,-,-(/•) = 9{r - aij) H -5{r - o^) , (5.1) 



'17 V 



where 9 is the Heaviside step function and 5 the Dirac delta function. 

When studying percolation problems in the continuum is useful to rewrite the Boltzmann 
factor as the sum of two terms [23113 e ij( r ) — e ij( r ) + e tj( r )i where 

4(r) = 9(r-a i3 ), (5.2) 
e+.(r) = ^±5(r-a t] ). (5.3) 

The corresponding Mayer functions will be fij(r) = f*Ar) + fij(r), with 

fi(r) = 4(r)-l, (5.4) 
f$(r) = e+(r). (5.5) 

The procedure to obtain equations of connectedness and blocking functions from the usual 
pair correlation functions and direct correlation functions is best described through the use of 
graphical language. If we substitute and f£ bonds for bonds in the density expansions 
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for these functions, then the connectedness functions, which we will indicate with a cross 
superscript, are expressed as the sums of those terms that have at least one f£ bond path 
connecting the two root vertexes. The sums of the remaining terms in the expansions give the 
blocking functions. 

The percolation threshold corresponds to the existence of an infinite cluster of particles and 
is given by the divergence of the mean cluster size fIG\ ITS] 



where hf-{r) is the pair connectedness function (related to the joint probability of finding a 
particle of species i and a particle of species j at a distance r and that these two particles are 
connected) and 



Since hfj(r) is related to the so called direct connectedness function c^(r) through an OZ 
equation, one can use Baxter formalism again, introducing a factor function qfj{r). If we now 

define Q +) ij(k) = 5y — 27r(p i pj) 1 / 2 q^(k), then it results that 



cluster 





(5.6) 




(5.7) 




(5.8) 



■in 



and thus 




(5.9) 



■in 



where 




(5.10) 



Clearly Q+\ m {0) diverges to infinity when |Q + (0)| = 0, and this relation defines the percolation 
threshold. 

Another interesting and related quantity is the average coordination number 




(5.11) 



V.l mMSA approximation 

The mMSA closure for cfJr) is 




(5.12) 
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On the other hand when r < we have e*Ar) = and f^{r) = e^r), so we must have 
exactly 



a 



yij{<Tij)S{r - Oij) r<a i:j . 



■i j 



Within the mMSA we have for the cavity function at contact [Zj 

Vij(crij) = 1 for all i, j . 



(5.13) 



(5.14) 



Following the same steps of Chiew and Glandt [HI HH] we then find (see Appendix IVHI for 
details) 



From which it follows 



Within Cases I, II, and III 



qtAr) = K ij d(r-L ij )d(a ij -r) 



Q+,ij(0) = 5y - 27r(p i p j ) 1/2 K ij cr j 



(5.15) 



(5.16) 



(5.17) 
(5.18) 
(5.19) 

Now from Eq. (|5.17|) follows that Q+.r,(0) is a 1-dyadic form. Using the properties of dyadic 
matrices (see Appendix IVII|) we then find 



Q+,y(0) = Sij + a+bf , 
aj = -lixp^fxiYi , 



IQ+(o)l 



at 1 + a+ ■ 



where 



|Q+(0)| = l + a+-b+ = 1-12£ 



1,2 



From Eq. ()5.10|) we find 



s m {0) 



IQ- 



x m (1 + a+ • b+) - 6+ J2 



(5.20) 



(5.21) 



(5.22) 



and from Eq. (|5.9j) 



24 ^1,1^0,1 , 144 6,2^,1 

^cluster J-~r* J - „ >. i" 



e 1 - i2fi, 2 e (i - 126,2) 



(5.23) 
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The percolation transition occurs when 

r <^> 3 

V 



a 3 ) ' M X M 2 



Wo 2 ) 



■V 



M. 



i] Case I , 
Case II , 
7] Case III 



(5.24) 



The threshold is independent of s at fixed rj for Case II, but lowers with increasing size 
polydispersity in Cases I and III. The curve is simply a straight line, as a consequence of the 
mean-field character of the mMSA (CO) closure. The qualitative result found with Cases I and 
III is however interesting. For the average coordination number we find from Eqs. (j5.1H) and 



iirpyixjXjKijcrij 



24 

r (a 3 ) 
^ {a){a 2 ) 
r {a 3 } 

At the percolation transition we then find 



(5.25) 



Case I 



Case ILIII 



2 Case I,III 
2/M 2 Case II . 



(5.26) 



Using Case IV Q + ij(0) turns out to be 3-dyadic; the percolation transition occurs when 
|Q+(0)|=0, i.e. 



1-3 

T 



s 2 {A + 7s 2 



o 



(5.27) 



■ (1 + 3s 2 + 2s 4 ) \tJ 16(1 + s 2 )(l + 2s 2 ) 2 \tJ 
The solution rj/r = p(s) such that p(0) = 1 is a monotonously decreasing function with 

lim p(s) = 0.756431 .... (5.28) 

Then with this Case we find that increasing the polydispersity the non-percolating region of 
the phase diagram diminishes. 

With Case V Q +i j(0) turns out to be 2-dyadic, and the percolation transition occurs when 



(a 3 ) i 2 



+ 



n 



M 2 V M X M 2 ) 2 ' 
which has the physical behavior already found with Cases I, II, and III. 



(5.29) 
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V.2 CI approximation with Case V 

As remarked, in Case V we can work out the percolation threshold equation even within the 
CI approximation. From Eq. (j5.13|) we have exactly 

hfjir) = —V^i^Kr ~ **) r < a tJ . (5.30) 

where (o^) is given by Eqs. IJ3.4|) . For the closure condition of the direct connectedness 
function we find again 

4 (r) = ft, (r) +fi(r)J2 ( r ) + % M E P^S™ W 

m m 

= r > aij , (5.31) 

since (r) = /y(r) = for r > a^. To determine gA(r) we then follow the same steps as for 
the mMSA case and we find 

g+(r) = Kfjy^ia^r - L ij )9(a ij - r) . (5.32) 

When we insert Kij from Eq. (|3.6p into the expression for Q +i j(0) [see Eq. (|5.16|) ] this 
becomes a 4-dyadic matrix whose determinant is 

6 

|Q+(0)| = 1 + ^^(3,77)/^, (5.33) 

i=l 

where the coefficients qi(s,rj) are given in Appendix IVII1 

The percolation threshold is the solution of |Q + (0)| = 0. This is an algebraic equation of 
order 6 in r. We can plot the correct root t(i]) for different values of polydispersity, as reported 
in Fig. El We see that increasing the polydispersity increases the non-percolating phase. One 
can clearly observe a clear improvement from the mMSA (CO) approximation although the 
7] — > limit is still qualitatively different from the PY one-component case. It would be 
interesting to study if the "true" percolation threshold passes through the origin (77 = 0, r = 0) 
(as occur in the CO or CI approximations) or has a finite limit (rj = 0, r = r ) (as it occur for 
monodisperse fluids in the PY approximation with tq = 1/12). Even if the Monte Carlo results 
of Ref. JOJE] are inconclusive in this respect, physically it is plausible to assume that at very 
low density the average number of bonds per particle is not sufficient to support large clusters 
at all and we would tend to favour the first scenario 3 . 

For the one-component system the average cluster size is 

S duster = 1 + ph + (0) 



l-pc + (0) [Q + (0)] 2 

(5.34) 



3 In this respect both CO and CI would be more precise than the PY closure and this is a remarkable feature. 
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The percolation transition occurs when rjy^ ci \a) 



= t or 



n 



2 (-3r 2 + v^r^Vl - 9r + 30r 2 ) 



(5.35) 



1 - 12r + 30r 2 



In Fig. |3]we compare our result for the one- component (s = 0) system with the PY result of 
Chiew and Glandt and the Monte Carlo simulation of Miller and Frenkel fUJ • 
The average coordination number becomes 



and at the percolation transition we find Z = 2. 

VI Phase equilibrium 

Phase equilibrium is another interesting aspect which can be analyzed in full details within 
our model. It was pointed out in Ref. jU] that the equation of state derived from the energy 
route for a one-component system of sticky hard spheres in the mMSA approximation is van 
der Waals like. The same holds true for the system studied with the CI approximation. It is 
worth stressing that the equation of state derived from the compressibility route cannot yield a 
van der Waals loop since from Eq. (j4.3J) [d{(3P) / dp]r > 4 . On the other hand the equation of 
state derived from the virial equation has been shown to diverge for the mMSA approximation 
[7 j and we anticipate that it also diverges for the CI approximation. This is the reason why we 
focus our analysis on the energy route in the present work. 

In this section we will find the binodal curves for the polydisperse system treated with the 
mMSA (CO) approximation and for the one-component system treated with the CI approxi- 
mation. The coexistence problem for a polydisperse system is, in general, a much harder task 
than its one-component counterpart, since it involves the solution of a large (infinite) number 
of integral non-linear equations. But we will see that since our excess free energy is expressed 
in terms of a finite number of moments of the size distribution function (a similar feature oc- 
curs for polydisperse van der Waals models [T2j, for polydisperse HS (THj and for Yukawa-like 
potentials ^ElES]) the coexistence problem can be simplified and becomes numerically solvable 
through a simple Newton- Raphson algorithm [see Eq. (|6.6p - (|6.8|) ]. The necessary formalism to 
this aim can be found in a recent review [IB] , and we will briefly recall it next. 

VI. 1 From a discrete to a continuous polydisperse mixture 

Consider a mixture made of p components labeled i = 1, . . . ,p, containing particles and 
with density p(°\ which separates, at a certain temperature r, into m daughter phases, where 
each phase, labeled a = 1, . . . ,m, has a number of particles and density p^ a \ Let the 
molar fraction of the particles of species % of phase a be , a = corresponding to the parent 
phase. At equilibrium the following set of constraints must be fulfilled: (i) volume conservation, 
(ii) conservation of the total number of particles of each species, (iii) equilibrium condition for 

4 Even though it may happen that one has loss of solution of Xia\ for certain values of the density, as 
occurs for the Percus Yevick closure [3]. 



Z 



(5.36) 
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the pressures 

p(«) (r, {x| a) }), and (iv) equilibrium condition for the chemical potentials 
Hi (t, p( a \ {x[ a '}). This set of constraints form a closed set of equations (see Appendix IVIII 
for details) for the (2 + p)m unknowns p(«), = N&/N(°\ and x\ a) with z = 1, . . . ,p and 
a = 1, ... ,m. Extension to the polydisperse case with an infinite number of components is 
achieved by switching from the discrete index variable i to the continuous variable a using the 
following replacement rule 

Xi -> F(cr)d<r , (6.1) 

where F(a)da is the fraction of particles with diameter in the interval (a, a + da). The function 
F(cr) will be called molar fraction density function or more simply size distribution function. 
Notice that, due to this replacement rule, we also have 

P (Q) (r,pM,K (a) }) -> P^{r,p^;[F^}) , (6.2) 

/iS a W a) .{*S o) }) - /i {Q W,p (a) ;[^ (a) ]) , (6-3) 

i.e. the thermodynamic quantities become functionals of the size distribution function and 
the equilibrium conditions (ii)-(iv) has to be satisfied for all values of the continuous variable 
a. The phase coexistence problem that now consists in solving the constraints (i)-(iv) for the 
unknowns p^ a \ x {a \ and F^ a \a) for a = 1, . . . , m, turns out to be a rather formidable task 
hardly solvable from a numerical point of view. Fortunately, as outlined in the next subsection, 
for our model a remarkable simplification occurs. 

VI. 2 Truncatable excess free energy 

As is described in the next subsection, the excess free energy of our system is truncatable: it 
is only a function of the three moments i = 1,2,3 of the size distribution function [see Eq. 
(I6~T21 for Case I, II, III, IV, and V treated with mMSA, and Eq. (l6~2T)j) for Case V treated 
with CI]. So we have the following simplification 

P H( T;P H.[ F (")]) ^ P (a) (r,^;{C ta) }). (6-4) 
^ a \tT,T,p^[F^}) - A* {a) (^r,pW;{e} o) }) , (6-5) 

where {^f^} is a short-hand notation for ^[ a \ ^\ £3 . It turns out that the two-phase (m = 2) 
coexistence problem, the one in which we are interested (we are thus concentrating on the 
high temperature portion of the phase diagram), reduces to the solution of the following eight 
equations in the eight unknowns p^\ p^ 2 \ {^f 1 ^}, and {^} 

= IpM jQ^\^p^\p^\p^{^},{^})F^\ayda 1 



z = l,2,3 a = 1,2, (6.6) 
1 = J Q ia Hcr } r } p^\p^\p^ ] {^\{^})F^(a)da , 

a = 1 or 2 , (6.7) 

i 5(1) (v (1) ;{d 1) }) = i 5(2, (v (2) ;{C (2) }), (6.8) 
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with 



p (a) Q (a) = p (0) 



( p (l) _ p (0)) + ( p (0) _ p (2)) e /3A^- 



(6.9) 



and 



A// 



z{2) 



(*,T, P WaF {2) ))-» eXCW (cT,T,pM-,[FW)) 



(6.10) 



where we indicate with the superscript exc the excess part (over the ideal) of the chemical 
potential. For a complete derivation of Eqs. ()6.6j) - (j6.8|) see Appendix IVIII 

VI. 3 Thermodynamic properties 



In order to obtain the equation of state of our model Eq. (j2.1)l from the energy route, one 
exploits the following exact result [I] 



d{(3A exc /N) 



Or 



gij(r)r 2 dr 



2-np XiXj I 

2np2_^XiXj / - eij{r)yij{r)r 2 dr 
2np J 10 L D Z ViA r ) r2dr ■ 



Upon taking the sticky limit we find 

d((3A exc /N) 
~dr 



Tij 



(6.11) 



VI. 3.1 mMSA approximation 

Within the mMSA approximation the partial cavity functions at contact are all equal to 1 so 
from Eq. (|6.11|) . after integration over r from r = oo (hard sphere case), we find 



exc a exc \ 

SHS S±HS) 



Case I , 
Case II, III , 



N 



—66 

Ti 

— t (366 + 66) Case IV 
r 4 

-?K 66+ I) CaseV - 



(6.12) 



The pressure can be found, from /3P/p = i]d(j3A/N)/drj 



^(3[P SHS (r,p; {&}) - P HS (r,p; {6})] = 0{A ^ 



6, 



(6.13) 
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where for Phs we use an equation due to Boublik, Mansoori, Carnahan, Starling, and Leland 
(BMCSL) fI7\ 12*5] which reduces to the Carnahan-Starling one when s = 0, 



7T 



-PPHs(T,p;{Zi}) 

6 



6 



6 



1-6 



+ 



+ 3- 
i 

?yrj 



66 



6$ 



1-7] (1-77) 2 M 2 



i - 6) 2 



+ 



;i-6) 3 a 



3 V 2 



1(1 - V f (l-vf 



6) 3 

Mi 
Mf 



(6.14) 



The excess free energy of the polydisperse hard sphere system is obtained integrating (Zhs~ 
l)/rj over rj, from 77 = 0, and recalling that the excess free energy is zero when rj = 0. We then 
find HI] 



f3A 



exc 
HS 



rj 



Mi 3r/ 1 



iV 



+ 



;i-r/) 2 M| l-r]M 2 



+ 



Mi 



M| 



66 



66(1-6 



6(1-6 



+ 



66 2 



ln(l — 77) 



1 ln(l - 6 



(6.15) 



Note that both A e ^ s and A e £ c s depend upon only a finite number of moments £ v , and this is 
the crucial feature for the feasibility of the phase equilibrium, as remarked. 
For the chemical potential = d((3A/V)/dpi we find after some algebra 



r, p; {6}) = (A + A„[°]) + (jj$ 3 + a + 



v 2 +[p [ § s + Ap® 



a 



(6.16) 



where 



[2] 



[3] 



-ln(l-6) 

36/(1-6) 



•i) 

6 3 



ln(l - 6) + 36/(1 -6)+ 3^ /(1 -6 



ln(l-6)+(6-| 



2|)/(i-6) 3 



6 

/(i - 6) + (366 - 1 



(6.17) 
(6.18) 

(6.19) 



/(i-6) 2 + 



(6.20) 



and 



( HJ 



T A 

IJL 

I r26 2 



Case I , 
Case II, III , 
Case IV , 

Case V , 



(6.21) 
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[i] 



13g 

v° 

136 
t 4 



3«? 



Case I , 
Case II, III 
Case IV , 
Case V , 



(6.22) 



' t 4 
r 2 



Case I , 
Case II, III , 

Case IV , 

Case V , 



(6.23) 



e o 

o 



Case I , 
Case II, III 



— — Case IV , 
t 4 

Case V , 



(6.24) 



It is noteworthy that if we retain in the expression (|6.14|) for Phs, only the first term, 
then our Case IV coincides with the van der Waals model of Bellier-Castella et. al. ^7] with 
n = 1, 1 = 0, upon identifying 4r with the temperature used by these authors. 



VI. 3. 2 CI approximation with Case V 

In analogy with what we have done before, we now consider the CI approximation for Case V. 
Using Eq. O into Eq. f[6~TT|) 



d{f3A°*</N) _ 12 g 
dr r 

Integrating from r = oo we find 

o( Aexc _ Aexc\ 

N 



(<7 3 > 



(a 3 ) j 



Vl [{of , (a) (a 2 ) 



2r \{a 3 ) (a 3 ) 



+ 



i r {*){<?) , W) 3 , 



r 2 V4 (a 3 ) 2 4 (a 3 ) 2 
1 ( 1 (a) 6 + 1 <a» 2 > 



r 3 V.72 (a 3 ) 2 24 (a 3 ) 2 



+ 



(6.25) 



(6.26) 



18 



For this case we limit ourselves to study the coexistence problem for the one-component 
system. In table we compare the critical parameters obtained through the energy route for 
the mMSA, CI, PY approximations and MC simulation, for the one-component system. 

Notice that, as already pointed out in Ref. [7j, a density expansion of y(a) within the 
PY approximation gives to zero-order the y{o) of the mMSA approximation and to first-order 
the y(a) of the CI approximation (as should be expected comparing the density expansions of 
the closures corresponding to these approximations). So at low densities Zshs from mMSA, 
CI, and PY should be comparable. From table H] we see that the true critical parameters are 
between the PY and the CI ones. 

In Fig. 0] we depict the binodal curve obtained from the CI approximation for the one- 
component system and we compare it with the PY binodal curve (obtained from the energy 
route) jl] and the one resulting from the MC simulation of Miller and Frenkel [H]. Remarkably, 
the gas-liquid coexistence curve predicted by CI lies closer to the MC data than the one 
predicted by PY on the gas branch and further on the liquid branch. 

VI. 4 Numerical results 

In this section we describe the numerical results obtained from the solution of Eqs. (jfi.fijl - fjfi.Sj) 
for the SHS in the mMSA, through a Newton-Raphson algorithm. 

We first determined the cloud and shadow curves by solving Eqs. ()6.6|) - (J6.8|) for the particu- 
lar case in which we set p^ = p^ so that F^\a) = F^°'(a). The cloud curve p c (r) is such that 
the solutions p^(r), p^(r) of the full coexistence problem given by Eqs. (|6.6|) - (j6.8j) . for a fixed 
p(°) (the coexistence or binodal curves), have the property that for a certain temperature To, 
P^( r o) — Pc( r o) — P > i- e the density of phase 1 ends on the cloud curve. The shadow curve is 
the set of points p s (t) in equilibrium with the corresponding cloud curve, i.e. p( 2 )(r ) = p s (r ), 
the density of phase 2 ends on the shadow curve. The interception between the cloud and the 
corresponding shadow curve gives the critical point (r cr ,p cr ): when = p cr the two solutions 
P^'(t), p( 2 \r) meet at the critical point. 

In order to find the cloud and shadow curves we choose as the parent distribution (a) a 
Schulz distribution with (cr)=l, and the initial conditions, for the Newton-Raphson algorithm, 
in the high temperature r* and low polydispersity region. Our starting conditions for the 
solution are 



for a = 1,2, where p oc and p oc are the coexistence densities at a temperature r* for the one 
component system. Once the cloud and shadow curves are determined we proceed to find the 
coexistence curves for a given mother density. 

In Fig. El we depict the cloud and shadow curves obtained with our Case I for two rep- 
resentative values of polydispersity, s = 0.1 and s = 0.3. For comparison the coexistence 




.(«)(! + #(1 + 2*2) , 



{a) (l + sl), 



(6.27) 
(6.28) 



(6.29) 



(6.30) 
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curve of the one component system (s = 0) is also reported. As polydispersity increases, the 
critical point moves to lower densities and lower temperatures (r cr ~ 0.094, p cr ~ 0.249 at 
s = 0, r cr ~ 0.093, p cr ~ 0.24 at s = 0.1, and r cr ~ 0.085, p cr ~ 0.197 at s = 0.3). Let 
us now fix s = 0.3, a value corresponding to a moderate polydispersity. Again in Fig. EJwe 
depict three coexistence curves upon changing the mother density p^ = 0.08, p^ = 0.25, and 
p<® = 0.197 ~p cr . 

All these curves closely resemble the corresponding ones obtained for the polydisperse van 
der Waals model ^1 , in agreement with previous results. In Fig. El we show how the two 
daughter distribution functions (at s = 0.3 and P {0) = Per) 

differ from the parent Schulz 

distribution (a process usually called fractionation), for two different values of temperature 
r = 0.084 and r = 0.078. 

Next we consider differences in the critical behavior with respect to changement in the 
Case. In Fig. we show the cloud and shadow curves obtained using Case I, IV, and V at 
s = 0.3. While for Case I and V the critical point is displaced at lower temperature and lower 
density respect to the monodisperse system, the critical point of Case IV is displaced at higher 
temperatures and lower density. The cloud curves of Case II and III have a low density branch 
that does not meet the high density one as soon as s > 0; moreover, the cloud curve does not 
meet the corresponding shadow curve, so there is no critical point. We are not aware of similar 
features in other polydisperse models, although this is of course to be expected in other cases 
as well. 

VII Conclusions 

In this work we have performed an extensive analysis of the phase diagram for Baxter SHS 
model in the presence of polydispersity. In spite of its simplicity, this model has been proven 
to be extremely useful in the theoretical characterization of sterically stabilized colloids. These 
systems are, however, affected by intrinsic polydispersity in some of their physical properties 
(size, species, etc) and hence the effect of polydispersity on the corresponding theoretical mod- 
els cannot be overlooked and is then a rather interesting point to address. As only formal 
manipulations [3] can be carried out for the multicomponent Baxter SHS model within the PY 
approximation, we have resorted to a simpler closure (mMSA) to which the PY closure reduces 
in the limit of zero density and that was recently shown [7, to reproduce rather precisely many 
of the interesting features of its PY counterpart. Our analysis has also been inspired by re- 
cent results by Miller and Frankel ^T] who showed that Baxter SHS model coupled with PY 
closure reproduced fairly well their MC data in the one-component case. We have studied the 
effect of polydispersity on phase stability boundaries, the percolation phase transition, and the 
gas-liquid phase transition. We have considered five different cases of polydispersity. This is 
because there is no general agreement on the way in which adhesion forces are depending on 
the size of particles. Case I and II had already been discussed in previous work by us |B] , Case 
III is a variant of Case I, whereas a case similar to Case IV had been employed by Tutschka 
and Kahl ^2|- Finally Case V has been specifically devised to cope with approximation CI. In 
spite of the apparent redundancy of all these sub-cases, we believe that each of these examples 
has a reasonable interest on its own, and hence we have included them all in our discussion. 
We studied the instability boundaries and the two-phase coexistence problem of polydisperse 
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SHS system in the mMSA (CO). The next level of approximation (CI) would still be feasible, 
but significantly more lengthly. We have laid down the necessary formalism in Sections 11111 
and lV1.3~2l and tested its effect on the one-component case, by comparing the results against 
the PY approximation and MC data. We derived the percolation threshold of the polydisperse 
system both within mMSA (CO) closure (for all five Cases) and in the CI approximation (using 
CaseV). 

We found that the effect of polydispersity on the stability and phase boundaries slightly 
depends upon the chosen Case, but there are general features shared by all of them: polydis- 
persity renders the mixture more stable with respect to concentration fluctuations (in the small 
density region, see Fig. |2J) with the exception of Case II for which the stability boundary is 
independent from the polydispersity; Eqs. (|5.24|) . (|5.27|) . and ()5.29|) (in the mMSA), and Eq. 
(j5.rSHj) (in the CI), describe its effect on the percolation threshold (see Figs. E]and|HJ). Polydis- 
persity increases the region of the phase diagram where we have a non-percolating phase, with 
the exception of Case IV, for which the opposite trend is observed, and of Case II for which the 
percolation threshold is independent from the polydispersity; polydispersity reduces the region 
of the phase diagram where we have a gas-liquid coexistence for Cases I and V, while the op- 
posite trend is observed for Case IV (see Fig. |7J). For Case II and III we obtained cloud curves 
with a gap at high temperature, moreover the cloud curve does not meet the corresponding 
shadow curve, so there is no critical point, as soon as polydispersity is introduced. 

In conclusion, the typical effect of polydispersity is to reduce the size of the unstable region, 
the percolating region, and the two-phase region of the phase diagram, although exceptions to 
this general rule have been observed for Case II, III, and IV. 

For the one-component case we also compared the percolation threshold and binodal curve 
obtained from the CI approximation with the results from the PY approximation fUll] and the 
results from the Monte Carlo simulation of Miller and Frenkel JT] (see Fig. 0]). The percolation 
threshold from CI appears to approach that from PY, but is still significantly different from the 
results from the Monte Carlo simulation (the zero density limit, on the other hand, appears to 
be more physically sound than the PY one, and this feature remains to be elucidated). The gas- 
liquid coexistence curve predicted by CI is better than the one given by PY on the gas branch 
and worse on the liquid branch. Table |l] shows how the true (from the Monte Carlo simulation 
of Miller and Frenkel JT]) critical temperature and density for the gas-liquid coexistence should 
lay between the ones predicted by PY and the ones predicted by CI. 

Future developments of the present work can be envisaged along the following lines: (i) as 
pointed out in |221 on defining = Yl m x m / 'S C c(fy and = Yl m x m /[(pk B TK T )S C c(0)], the 
condition ipQ > is necessary but not sufficient for the material stability of the system and 
the condition ip^ > is necessary but not sufficient for the mixed material and mechanical 
stability. It could happen that those two conditions are satisfied but the system is nonetheless 
unstable as occur for example in the binary mixture studied by Chen and Forstmann ■ Even 
though a characterization of the instability boundary in the spirit of Chen and Forstmann seems 
unattainable for a polydisperse system, it would be desirable, in the future, a more precise 
location of the instability boundaries. Moreover the way we found the instability boundaries 
for the polydisperse system was to start from the instability condition valid for a discrete 
mixture and take the limit of a continuous mixture on the instability boundaries of the discrete 
mixture. It would be interesting to compare our analysis with the one given by Bellier-Castella 
et. al (see section II C in Ref. [T7] ) who take the continuous limit from the outset; (ii) 
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all the percolation thresholds that we have calculated have a low density branch that enters 
the gas-liquid coexistence region. The same feature is observed for the one-component system 
studied through Monte Carlo simulation [TlU lllj. While it is clear that continuum percolation is, 
strictly speaking, not a thermodynamic phase transition, one could expect, from a "dynamical" 
point of view, an interference between the formation of infinite clusters of particles and phase 
separation, and a clarification of this point would have interesting experimental applications; 
(iii) the polydisperse system is expected to display, in the low temperature region, other critical 
points signaling the onset of m > 2 phase coexistence [TJ] , and it would be interesting to study 
their evolution with polydispersity for our system. 



Appendix A: Determinant and inverse of a dyadic matrix 

Given the n-dyadic matrix of Eq. ()4.6|) . its determinant is 



Q 



1 + A« 


B (D 


Ad) 


B (2) .. 


Ad) 


■ Bd) 


A (2) 


B d) 


1 + A( 2 ) 


• B( 2 ) •■ 


A (2) 


■ Bd) 


aw 




A (n) 


. B ( 2 ) .. 


■ 1 + A<") 


Bd) 



(A.l) 



Furthermore, any dyadic matrix Q admits analytic inverse for any number p of components, 
with elements given by 
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B: Derivation of Eq. (5.15) 







(A.2) 



The closure condition (j5.12J) justify the usual generalized Wiener-Hopf factorization [3*T] 

/oo 
dtq+Mq+Zir + t) , (B.3) 
,,, -mi 

/OO 
dtqL(t)(r - t)h + m] {\r - t\) , (B.4) 
,,, -im 

(B.5) 

where r > L^-, the prime denotes differentiation, and qfj{r) are real functions with support on 
[Lij,<Tij] and zero everywhere else. 

Let us determine qfj{r). Using the exact condition ()5.13)1 in Eq. ()B.4)) we find for < r < 



q+'(r) = -K ij 5{\r\-a ij ) + 2^ t 



K. 



dtqUt)(r-t)-^5(\r-t\-a mj ). (B.6) 



a 



rnj 
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The second term on the right end side is equal to 2tt Ylm Pm,citm{ r ~ a mj)K m j which is zero when 
r < <7jj. So we simply have 

Integrating this equation gives Eq. (|5.15j) . 



Appendix C: Coefficients of Eq. ( 15.331) 

The coefficients in Eq. fl5.33j) are as follows 

, , y (2 + 577) (l + 3s 2 + 23 4 ) 3 

qi{s,r]) = 7 , (C8) 

V ' 2(l + s 2 ) 3 (l + 2s 2 ) 4 



"■ ■ 3 2 ) 3 (l + 2s 2 ) 4 

, , '/ 1 ~~ ~ r ['■> '/ 1 1 ~ r 7j ~~ 5] s 2 — 2 s 4 } , 

= ~ — ^77— 2 . 3 > ( C - 10 ) 

3 ) (1 + 2 s z ) 

, , '/ - I- - -^/ + (4 + 777) s 2 } 

*M = ^ , ^277—^4 , (Cll; 



if {-4 + [77 (2 + 77) — 5] s 2 } (1 + 3 s 2 + 2 s 4 ) z 

4(l + s 2 ) 3 (l + 2s 2 

r? 2 {-2 + [6 77 (1 + 77) - 5] s 2 - 2 s 4 } 

24 (1 + s 2 ) (l + 2s 2 

77 3 s 2 [2 + 577 + (4 + 7r]) s 2] 

96(l + s 2 ) 2 (l + 2s 2 

q 5 (s,v) = 0, (C.12) 

4 4 

g 6 (s,77) = ^-4 j. (C.13) 

V " 2304(l + s 2 ) 3 (l + 2s 2 ) 4 V ; 

Appendix D: Phase coexistence conditions 

In this Appendix we give a complete derivation of Eqs. (|fj.fjj) - (|fj.8jl in the main text. 

Consider a p— component mixture. Each species i has number density p\ = N- ^ /V^°\ 
where N- ^ is the number of particles of type i and the volume of the system. 

We assume that at a certain temperature r the system separates into m daughter phases, 
where each phase a = 1, . . . , m is characterized by a volume and a number of particles of 
species 1, iv> . 

At equilibrium the following set of constraints must be fulfilled: 
(1) volume conservation 



= j2 y {a) , (D.i4) 



(2) conservation of the total number of particles of each species 

m 

NP=Y^Ni a \ i = l,...,p, (D.15) 
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(3) equilibrium condition for the pressures 

P {a \r, V {a \ {N^}) = pW(r, V {/3 \ {N^}) , (D.16) 

(4) equilibrium condition for the chemical potentials 

^\rM a \{N^})=^\rM"\{Nf ] }) , i = l,... lPl (D.17) 

where {N^} is a short-hand notation for iVf, . . . , N£. 

It is convenient to use the following set of variables: r; p^ = JV^/V^; scj"' = iVj /N^ a \ 

i = l,...,p with M Q ) = J2i N i a) - Introducing = N^/N® Eqs. (jDTT4|l - (|DT7jl can be 
rewritten as follows 



1 




(D.18) 


r (0) 


= E^ (a) > 


(D.19) 


i a) }) 


= P^(r,p^,{xf}), 


(D.20) 


i a) }) 


= pfW^f}), 


(D.21) 



/4 Q W Q) , 

with the normalization condition 

J> 4 (Q) = 1, a = l,...,m. (D.22) 

Eqs. flD.18|) - ()D.22|) form a set of closed equations for the (2 + p)m unknowns 
with i = 1, . . . ,p and a = 1, . . . , m. Notice that when m = p + 1 the densities of the various 
phases p^ will be independent of p(°\ since relations (jD.20|) . ()D.21|) . and ()D.22|) form a closed 
set of equations for the unknowns 

In the continuous polydisperse limit (p — > oo) we have to take into account the substitution 
rule (jbM|) . Then the thermodynamic quantities will be rewritten as in Eqs. (J6.2j) and (|6.3|) . and 
Eqs. dnHt - dESl) become 



E^ (Q) ' (D.23) 
F (0) (a) = , (D.24) 



p(o) ~ p {a) ■ 



pW(r,/)W; [F^]) = P^\r,p^; [F^]) , (D.25) 
^\a,T,pW;[F^)) = ^\a,r,p^;[F^)), (D.26) 



with the normalization condition 
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Integrating Eq. (|D.24|) over a and using Eq. (|D.27|) we find 

£V a) = l. (D.28) 

a 

The set of Eqs. ()D.23j) - (jD.27j) form a closed set of equations for the unknowns and 
F^(cr) with a = l,...,m. Notice that, due to the substitution rule (jfi.lj) . sum over i be- 
comes integration over a and thermodynamic quantities become functionals of the distribution 
function. We have indicated such dependence with square brackets. 



Two-phase coexistence 

Let us now specialize ourselves to the two-phase (m = 2) coexistence. We are thus concentrating 
on the high temperature portion of the phase diagram, since coexistence with m > 2 (Gibbs 
phase rule does not restrict the value of m in a system of infinitely many species) is expected 
to occur at low temperatures. From Eqs. (jD.28|) and (jD.23|) we find 

(1) P (0) - P (2) P (1) 

X = p (l) _ p (2) p (0) ' ( D ' 29 ) 

(2) - P (0) P (2) 

X2 = p (l) _ p {2) ^0) ■ (D ' 30) 

Notice that x^- 1 ' and x^ must be positive. So if p^ 1 ' < p^ 2 \ then p^ must lie between p^ 
and p( 2 \ if p^ < p^\ it must lie between p^ and p^\ Substituting these expressions in Eq. 
(lD~2lj) we find 

p(2)p(2 , = P <». F <».^-^ + /"F".^-^ . (D.31) 

Next we divide the chemical potentials in their ideal and excess components /i = p ld + p exc 
where 

(3p id(a \a } r } p^; [FW]) = lnfAV^^^H^)] , (D.32) 
with A being the de Broglie thermal wavelength. Now Eq. (|D.26j) becomes 

= F {2 \a)^e^\ (D.33) 
pW 

A/i e * c = /i exc(2) (a,r,p( 2 );[F( 2 )])-/i exc(1) (a,r,p«;[F«]) . (D.34) 
From Eqs. (JD~3TJ) and (fTT33j) we find 

F^\a) = F^(a)Q^(a,r,p^,p^,p^; [F< 2 >]) , (D.35) 

where the are defined by Eq. (|6.9jl . 

Formally the set of Eqs. (JEHU), (lD~33j) . (lD~25|) with a = l,/3 = 2, and (jEZTj) with a = 1 
or 2, form a closed set of equations for the unknowns p^ l \p^ 2 \ F^ l '[a) and F^ 2 \a). In our case 
the free energy of the system [Case I, II, III, IV, and V treated with mMSA, see Eq. (JbM2|) . 
or Case V treated with CI, see Eq. (|6.26|) ] is truncatable: it is only a function of the three 
moments £j, i = 1,2,3 of the size distribution function F. So the problem is mapped in the 
solution of the 8 Eqs. in the 8 unknowns p«, , ff 5 , &\ , ^ 2) , ^ 2) , ^ . 
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LIST OF FIGURES 



Fig. Schematic diagram showing the area of the contact surface between a particle of species 
i and a particle of species j. 

Fig. 121 Curves for the onset of phase instability (the fluid is stable above the curves shown) as 
obtained from the mMSA approximation for a monodisperse s = system, and for a 
polydisperse system with s = 0.2, and polydispersity chosen as in Cases I, II, and III 
[see Eq. (j4.10|) ]. We also show for the one-component system the curve for the onset of 
mechanical instability predicted by the CI approximation [see Eq. (j4.12|) ] and the one 
predicted by the PY approximation [see Eq. (j4.15j) . 

Fig. |3] Dependence of the percolation threshold, as calculated from the CI approximation using 
Case V (see section E2J), from the polydispersity. 

Fig. 0] Binodal curve and percolation threshold [see Eq. (|5.35[) ]. for a one-component system, 
in the CI approximation. For comparison we also show the percolation threshold of the 
Percus-Yevick approximation (which exists for r > 1/12), the one from the Monte 
Carlo simulation of Miller and Frenkel [H] (circles are the simulation results and the fit, 
the dot-dashed line, is only valid for r > 0.095), the binodal curve of the Percus-Yevick 
approximation (from the energy route) and the binodal curve from the Monte Carlo 
simulation of Miller and Frenkel ^1] (points with errorbars are the simulation results and 
the fit, the dot-dashed line, is merely to guide the eye). 

Fig. El Cloud and shadow curves for Case I in the mMSA at two values of polydispersity: s = 0.1 
and s = 0.3. For the case s = 0.3 we also show three coexistence curves (continuous lines) 
obtained setting p^ = 0.08, p^ = 0.25, and p^ = 0.197 ~ p cr . For comparison the 
binodal of the monodisperse (s = 0) system has also been included. 

Fig. |B] Evolution of the size distribution of the coexisting phases F^(a) and F^(a), with tem- 
perature along the critical binodal of Fig. EJ(s = 0.3, p^ = 0.197 ~ p cr ). For comparison 
also the parent Schulz distribution is shown (continuous line). 

Fig. Cloud and shadow curves for Case I, IV, and V in the mMSA at s = 0.3. For comparison 
the binodal of the monodisperse (s = 0) system has also been included (continuous line). 
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LIST OF TABLES 

Table HI For the one-component system, we compare the critical parameters obtained from the 
mMSA, CI, and PY [I] approximations with the ones from the Monte Carlo simulation 
of Miller and Frenkel ^1] . 
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Tc Vc (Zshs)c 

mMSA 0.0943 0.13 0.36 

CI 0.1043 0.14 0.37 

PY 0.1185 0.32 0.32 

MC 0.1133 0.27 



Table I: R. Fantoni, D. Gazzillo, and A. Giacometti 
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a. a. 




Figure 1: R. Fantoni, D. Gazzillo, and A. Giacometti 
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Figure 2: R. Fantoni, D. Gazzillo, and A. Giacometti 
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Figure 5: R. Fantoni, D. Gazzillo, and A. Giacometti 
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Figure 7: R. Fantoni, D. Gazzillo, and A. Giacometti 
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